The aims of this notebook are to (1) download and read the 10X Multiome data for two mantle cell lymphomas (MCL) presented in tonsils, (2) perform basic quality control to exclude cells and doublets, (3) reduce the dimensionality of the data (PCA, UMAP), (4) cluster cells into meaningful cell states, and (4) annotate these cell states based on their top marker genes.
library(Seurat)
library(tidyverse)
library(scDblFinder)
library(DT)
library(ggpubr)
library(here)
library(glue)
library(Polychrome)
plot_histogram_qc <- function(df, x, x_lab) {
df %>%
ggplot(aes_string(x)) +
geom_histogram(bins = 100) +
labs(x = x_lab, y = "Number of Cells") +
theme_pubr()
}
process_seurat <- function(seurat_obj, dims = 1:30, reduction = "pca") {
seurat_obj <- seurat_obj %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
RunUMAP(dims = dims, reduction = reduction) %>%
FindNeighbors(dims = dims, reduction = reduction)
seurat_obj
}
# Thresholds
min_lib_size <- 900
min_n_genes <- 700
max_pct_mt_102 <- 27.5
max_pct_mt_413 <- 20
chrY_cutoff_102 <- -0.32
chrY_cutoff_413 <- -0.1
The expression and accessibility matrices for these two MCL patients are available in this Zenodo respository. We can download the data from the terminal using wget as follows:
cd MCL
mkdir -p data
wget -P data https://zenodo.org/record/6678331/files/TonsilAtlasCellRangerOuts.zip
cd data
unzip TonsilAtlasCellRangerOuts.zip
Then, we can proceed to load the data into R. Note that in this project we only focus on the RNA matrices. The accessibility information will be used in another project. We will also proceed to read the metadata:
# Read metadata
metadata <- read_csv(
here("MCL/data/TonsilAtlasCellRangerOuts/MCL/metadata/cellranger_arc_metadata.csv"),
col_names = TRUE
)
metadata <- metadata[str_detect(metadata$library_name, "pos"), ]
metadata <- metadata[metadata$type == "RNA", ]
DT::datatable(metadata)
gem_ids <- unique(metadata$gem_id)
# Read expression matrices and convert them to Seurat objects
path_to_mats <- here(glue("MCL/data/TonsilAtlasCellRangerOuts/MCL/BCLLATLAS_64_65/{gem_ids}/filtered_feature_bc_matrix"))
names(path_to_mats) <- gem_ids
se_l <- map(gem_ids, \(x) {
counts <- Read10X(path_to_mats[x])$`Gene Expression`
seurat_obj <- CreateSeuratObject(counts)
seurat_obj$gem_id <- x
new_metadata <- seurat_obj@meta.data %>%
rownames_to_column("barcode") %>%
left_join(metadata, by = "gem_id") %>%
as.data.frame()
rownames(new_metadata) <- new_metadata$barcode
seurat_obj@meta.data <- new_metadata
seurat_obj
})
metadata <- as.data.frame(metadata)
rownames(metadata) <- metadata$gem_id
names(se_l) <- metadata[gem_ids, "donor_id"]
se_l <- map(se_l, \(seurat_obj) {
seurat_obj$pct_mt <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")
seurat_obj$pct_ribo <- PercentageFeatureSet(seurat_obj, pattern = "^(RPS|RPL)")
seurat_obj
})
# QC metrics across samples
qc_df <- map_df(se_l, \(seurat_obj) seurat_obj@meta.data)
qc_metrics_rna <- c("nCount_RNA", "nFeature_RNA", "pct_mt")
qc_ggs_rna <- map(qc_metrics_rna, \(x) {
p <- ggplot(qc_df, aes_string("library_name", x, fill = "donor_id")) +
geom_violin() +
xlab("") +
theme_classic() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "none"
)
p
})
qc_ggs_rna[[1]] <- qc_ggs_rna[[1]] +
scale_y_log10() +
ylab("Library Size (total number of UMI)")
qc_ggs_rna[[2]] <- qc_ggs_rna[[2]] +
ylab("Number of detected features")
qc_ggs_rna[[3]] <- qc_ggs_rna[[3]] +
ylab("Percentage of mitochondrial expression (%)")
qc_ggs_rna
## [[1]]
##
## [[2]]
##
## [[3]]
We can clearly see that both samples have a different distribution of QC metrics. Following the best practices, we will determine the thresholds for both samples separately. In particular, we see that sample M102 has a higher mitochondrial expression that one would expect for a single-nuclei assay. According to this technical note from Cellranger: Even in 10x assays with single nuclei as input (for example, Multiome assay), an elevated level of mitochondrial genes would mean that mitochondrial RNA could remain “stuck” to nuclear membranes or may get partitioned into GEMs and be detected by the gene expression biochemistry.
# Library size
(lib_size_hists <- map(se_l, \(seurat_obj) {
p <- plot_histogram_qc(
seurat_obj@meta.data,
x = "nCount_RNA",
x_lab = "Library Size (total number of UMI)"
)
p +
scale_x_log10() +
geom_vline(xintercept = min_lib_size, color = "red", linetype = "dashed")
}))
## $M413
##
## $M102
# Number of detected genes
(n_genes_hists <- map(se_l, \(seurat_obj) {
p <- plot_histogram_qc(
seurat_obj@meta.data,
x = "nFeature_RNA",
x_lab = "Number of detected features"
)
p +
geom_vline(xintercept = min_n_genes, color = "red", linetype = "dashed")
}))
## $M413
##
## $M102
# Percentage of mitochondrial expresion
pct_mt_hists <- map(se_l, \(seurat_obj) {
p <- plot_histogram_qc(
seurat_obj@meta.data,
x = "pct_mt",
x_lab = "Percentage of mitochondrial expression (%)"
)
p
})
pct_mt_hists[[1]] <- pct_mt_hists[[1]] +
geom_vline(xintercept = max_pct_mt_413, color = "red", linetype = "dashed")
pct_mt_hists[[2]] <- pct_mt_hists[[2]] +
geom_vline(xintercept = max_pct_mt_102, color = "red", linetype = "dashed")
pct_mt_hists
## $M413
##
## $M102
Let us visualize the QC metrics jointly:
# Number of genes vs total UMI colored by pct_mt
(joint_gg1 <- map(se_l, \(seurat_obj) {
p <- seurat_obj@meta.data %>%
ggplot(aes(nCount_RNA, nFeature_RNA, color = pct_mt)) +
geom_point() +
geom_vline(xintercept = min_lib_size, color = "red", linetype = "dashed") +
geom_hline(yintercept = min_n_genes, color = "red", linetype = "dashed") +
labs(
x = "Library Size (total number of UMI)",
y = "Number of detected features",
color = "Percentage of mitochondrial expression (%)"
) +
scale_color_viridis_c(option = "magma") +
theme_pubr()
p
}))
## $M413
##
## $M102
# % mitochondrial expression vs number of detected genes
joint_gg2 <- map(se_l, \(seurat_obj) {
p <- seurat_obj@meta.data %>%
ggplot(aes(nFeature_RNA, pct_mt)) +
geom_point(size = 0.25, alpha = 0.5) +
geom_vline(xintercept = min_n_genes, color = "red", linetype = "dashed") +
labs(
x = "Number of detected features",
y = "Percentage of mitochondrial expression (%)"
) +
theme_pubr()
p
})
joint_gg2$M413 <- joint_gg2$M413 +
geom_hline(yintercept = max_pct_mt_413, color = "red", linetype = "dashed")
joint_gg2$M102 <- joint_gg2$M102 +
geom_hline(yintercept = max_pct_mt_102, color = "red", linetype = "dashed")
joint_gg2
## $M413
##
## $M102
# 102
is_low_quality_102 <-
se_l$M102$nCount_RNA < min_lib_size |
se_l$M102$nFeature_RNA < min_n_genes |
se_l$M102$pct_mt > max_pct_mt_102
table(is_low_quality_102)
## is_low_quality_102
## FALSE TRUE
## 10103 965
se_l$M102 <- se_l$M102[, !is_low_quality_102]
# 413
is_low_quality_413 <-
se_l$M413$nCount_RNA < min_lib_size |
se_l$M413$nFeature_RNA < min_n_genes |
se_l$M413$pct_mt > max_pct_mt_413
table(is_low_quality_413)
## is_low_quality_413
## FALSE TRUE
## 11074 2001
se_l$M413 <- se_l$M413[, !is_low_quality_413]
We follow the dimensionality reduction pipeline from Seurat. This includes (1) finding highly variable genes, (2) scaling data, (3) running principal component analysis (PCA), (4) running UMAP based on PC coordinates, (5) embedding cells in a K-nearest neighbors graph (representing cells as principal components).
se_l <- map(se_l, process_seurat)
map(se_l, DimPlot, reduction = "umap")
## $M413
##
## $M102
We follow the “Doublet detection by simulation” section from the book “Orchestrating Single-cell Analysis With Bioconductor”:
se_l <- map(se_l, \(seurat_obj) {
sce <- as.SingleCellExperiment(seurat_obj)
seurat_obj$doublet_score <- computeDoubletDensity(sce, dims = 30, k = 20)
seurat_obj
})
map(se_l, \(seurat_obj) {
p <- FeaturePlot(seurat_obj, "doublet_score") +
scale_color_viridis_c(option = "magma")
p
})
## $M413
##
## $M102
The first step will be to cluster cells and identify clusters of doublets using multiple sources of evidence (QC metrics, markers, scDblFinder). Second, we will cluster and annotate cells to a specific cell state.
se_l$M102 <- FindClusters(se_l$M102, resolution = 0.4)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 10103
## Number of edges: 360948
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8462
## Number of communities: 9
## Elapsed time: 1 seconds
cols <- kelly.colors(length(levels(se_l$M102$seurat_clusters)) + 1)
names(cols) <- NULL
cols <- cols[2:length(cols)]
DimPlot(se_l$M102, cols = cols)
We see that cluster 1 seems to correspond to the cluster that also has a high doublet score. Let us assess the distribution of QC metrics:
vln_qc_clusters <- VlnPlot(
se_l$M102,
c("nCount_RNA", "nFeature_RNA", "pct_mt", "doublet_score"),
pt.size = 0,
ncol = 2,
cols = cols
)
vln_qc_clusters[[1]] <- vln_qc_clusters[[1]] +
scale_y_log10()
vln_qc_clusters
In fact, we see that both cluster 1 and cluster 6 have a disproportionately high number of counts, detected genes, and doublet score. We proceed to find markers for these clusters:
markers_1_6 <- map(c("1", "6"), \(x) {
df <- FindMarkers(
se_l$M102,
ident.1 = x,
only.pos = TRUE
)
head(df, 40)
})
markers_1_6
## [[1]]
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## EEPD1 4.119048e-102 0.3016145 0.472 0.180 1.507613e-97
## SLC44A5 6.159026e-94 0.3499780 0.317 0.105 2.254265e-89
## GFOD1 2.689827e-90 0.2545392 0.471 0.187 9.845035e-86
## ATRNL1 3.621253e-80 0.2787145 0.475 0.202 1.325415e-75
## AC090125.1 8.104562e-77 0.2591640 0.401 0.163 2.966351e-72
## ZNF804A 1.072725e-76 0.4079666 0.999 0.944 3.926282e-72
## LY75 1.863314e-71 0.2524651 0.580 0.272 6.819915e-67
## PDE3B 1.295099e-70 0.2598599 0.691 0.361 4.740191e-66
## EXT1 1.777397e-68 0.4033100 0.915 0.645 6.505451e-64
## ADARB2 1.597433e-65 0.3141427 0.670 0.365 5.846764e-61
## PLCL2 2.586330e-65 0.3983377 0.980 0.807 9.466227e-61
## IGF2BP3 5.507840e-64 0.2675866 0.755 0.419 2.015925e-59
## ROR1 6.122316e-64 0.3488344 0.993 0.914 2.240829e-59
## HS2ST1 8.852155e-63 0.2678228 0.658 0.349 3.239977e-58
## TMEM65 1.554270e-62 0.2811789 0.707 0.384 5.688783e-58
## LTBP1 5.804061e-61 0.2826305 0.710 0.389 2.124344e-56
## MAN1A1 2.896576e-58 0.3578664 0.944 0.705 1.060176e-53
## MSI2 8.199324e-57 0.3115818 0.994 0.918 3.001034e-52
## OSBPL3 1.967148e-56 0.3376223 0.922 0.657 7.199958e-52
## AL136962.1 2.220951e-55 0.3249047 0.834 0.540 8.128903e-51
## KHDRBS2 1.146607e-54 0.2590180 0.743 0.439 4.196695e-50
## PSD3 2.425444e-54 0.2600765 0.786 0.458 8.877367e-50
## PLCB1 4.209754e-54 0.3420264 0.710 0.424 1.540812e-49
## CMSS1 1.376715e-53 0.2512005 0.766 0.466 5.038913e-49
## CCDC88A 1.049070e-52 0.3145087 0.812 0.498 3.839702e-48
## GALNT1 1.288633e-49 0.2683721 0.792 0.466 4.716525e-45
## TCF4 1.347879e-49 0.3171343 0.968 0.807 4.933372e-45
## JAZF1 5.592023e-49 0.3383225 0.961 0.775 2.046736e-44
## GAB2 2.950968e-47 0.2881779 0.988 0.867 1.080084e-42
## MALT1 1.671822e-46 0.2516282 0.884 0.592 6.119037e-42
## RGS12 1.996903e-46 0.2500715 0.904 0.638 7.308864e-42
## ITFG1 5.612126e-46 0.3003210 0.864 0.595 2.054094e-41
## ZNRF2 2.118019e-44 0.2521834 0.827 0.521 7.752162e-40
## TSHZ2 4.747068e-43 0.2744857 0.950 0.710 1.737474e-38
## TRPS1 6.198662e-43 0.2680878 0.927 0.694 2.268772e-38
## ITPR2 5.034643e-41 0.2743261 0.934 0.693 1.842730e-36
## IMMP2L 5.191785e-41 0.2521681 0.899 0.654 1.900245e-36
## FNDC3A 6.565079e-41 0.2582182 0.900 0.623 2.402885e-36
## TMEM163 5.364602e-40 0.2707752 0.956 0.765 1.963498e-35
## MCTP2 1.108299e-39 0.2744925 0.872 0.595 4.056486e-35
##
## [[2]]
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## CDCA8 0.000000e+00 0.4979704 0.289 0.006 0.000000e+00
## KIF2C 0.000000e+00 0.7385616 0.403 0.014 0.000000e+00
## NUF2 0.000000e+00 0.9097619 0.505 0.025 0.000000e+00
## ASPM 0.000000e+00 1.7927932 0.669 0.016 0.000000e+00
## KIF14 0.000000e+00 0.9325764 0.462 0.014 0.000000e+00
## CENPF 0.000000e+00 2.0241955 0.695 0.041 0.000000e+00
## RRM2 0.000000e+00 1.5572361 0.649 0.018 0.000000e+00
## NCAPH 0.000000e+00 0.6819171 0.364 0.009 0.000000e+00
## CKAP2L 0.000000e+00 0.6660459 0.321 0.003 0.000000e+00
## SPC25 0.000000e+00 0.5726783 0.275 0.004 0.000000e+00
## KIF15 0.000000e+00 0.9093094 0.446 0.018 0.000000e+00
## POLQ 0.000000e+00 1.6414247 0.695 0.049 0.000000e+00
## ECT2 0.000000e+00 0.8235383 0.361 0.011 0.000000e+00
## NCAPG 0.000000e+00 0.8516475 0.446 0.011 0.000000e+00
## HIST1H3B 0.000000e+00 1.0575366 0.426 0.009 0.000000e+00
## HIST1H1B 0.000000e+00 1.5948073 0.636 0.052 0.000000e+00
## KIFC1 0.000000e+00 0.7923254 0.380 0.008 0.000000e+00
## ANLN 0.000000e+00 0.7054958 0.348 0.006 0.000000e+00
## NCAPG2 0.000000e+00 1.5227158 0.708 0.058 0.000000e+00
## ESCO2 0.000000e+00 0.5126381 0.243 0.004 0.000000e+00
## MELK 0.000000e+00 0.8686628 0.439 0.014 0.000000e+00
## CDK1 0.000000e+00 0.5857226 0.298 0.004 0.000000e+00
## KIF11 0.000000e+00 1.1009219 0.538 0.019 0.000000e+00
## MKI67 0.000000e+00 2.0348210 0.823 0.016 0.000000e+00
## E2F8 0.000000e+00 0.5736032 0.249 0.002 0.000000e+00
## CDCA5 0.000000e+00 0.6658284 0.331 0.011 0.000000e+00
## DIAPH3 0.000000e+00 1.9461695 0.767 0.041 0.000000e+00
## DLGAP5 0.000000e+00 0.4790188 0.246 0.003 0.000000e+00
## BUB1B 0.000000e+00 0.7597472 0.384 0.011 0.000000e+00
## NUSAP1 0.000000e+00 1.5931659 0.757 0.037 0.000000e+00
## KIF23 0.000000e+00 0.6590418 0.305 0.004 0.000000e+00
## SHCBP1 0.000000e+00 1.0745399 0.541 0.017 0.000000e+00
## AURKB 0.000000e+00 0.5470624 0.272 0.003 0.000000e+00
## TOP2A 0.000000e+00 0.9309511 0.413 0.011 0.000000e+00
## KIF18B 0.000000e+00 0.5943851 0.275 0.004 0.000000e+00
## BIRC5 0.000000e+00 0.8907318 0.430 0.005 0.000000e+00
## TPX2 0.000000e+00 0.7740614 0.390 0.011 0.000000e+00
## GTSE1 0.000000e+00 0.9687339 0.479 0.006 0.000000e+00
## HJURP 6.782255e-301 0.3503057 0.187 0.002 2.482373e-296
## ESPL1 1.176614e-297 0.5432808 0.282 0.008 4.306524e-293
Cluster 1 does not have any specific marker (all markers have a very low log2 FC). Cluster 6, on the other hand, consists of proliferative cells, which is consistent with the high number of counts and features. Let us calculate a per-cell S phase and G2M phase score:
se_l <- map(se_l, \(seurat_obj) {
seurat_obj <- CellCycleScoring(
seurat_obj,
s.features = cc.genes.updated.2019$s.genes,
g2m.features = cc.genes.updated.2019$g2m.genes
)
seurat_obj
})
map(se_l, FeaturePlot, features = c("S.Score", "G2M.Score"))
## $M413
##
## $M102
Finally, we run the findDoubletClusters from scDblFinder, which identifies clusters with an expression profile lying between two other clusters:
sce_M102 <- as.SingleCellExperiment(se_l$M102)
dbl_out <- findDoubletClusters(
sce_M102,
clusters = sce_M102$seurat_clusters
)
dbl_out
## DataFrame with 9 rows and 9 columns
## source1 source2 num.de median.de best p.value lib.size1 lib.size2 prop
## <character> <character> <integer> <numeric> <character> <numeric> <numeric> <numeric> <numeric>
## 7 8 5 26 1530.0 ATXN1 2.53992e-09 1.049385 1.422903 0.00772048
## 0 7 1 86 535.5 MT-CO2 5.23716e-22 0.654574 2.219932 0.40473127
## 2 7 1 86 745.5 LINC01250 1.68773e-40 0.652915 2.214306 0.14668910
## 3 6 5 114 609.0 TSHZ2 3.06402e-35 1.844963 0.895737 0.09650599
## 5 7 3 171 936.0 LINC01250 7.21497e-62 0.702789 1.116399 0.07156290
## 8 7 4 236 2159.5 FYN 3.15284e-24 0.952939 1.681463 0.00752252
## 4 7 6 347 1112.5 ABTB2 1.33414e-86 0.566732 1.660963 0.08383648
## 1 7 6 395 1277.0 TSHZ2 7.55530e-32 0.294862 0.864174 0.15124221
## 6 8 1 625 1212.5 MKI67 2.45876e-75 0.358057 1.157174 0.03018905
We see that cluster 8 consists of T-cells (contamination / doublets):
FeaturePlot(se_l$M102, "CD3D")
VlnPlot(se_l$M102, "CD3D", cols = cols) + NoLegend()
markers_8 <- FindMarkers(
se_l$M102,
ident.1 = "8",
only.pos = TRUE,
logfc.threshold = 1
)
head(markers_8, 20)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## TGFBR3 0 2.023740 0.382 0.005 0
## CD2 0 1.824516 0.421 0.000 0
## CD247 0 2.866330 0.645 0.003 0
## STAT4 0 3.486855 0.724 0.017 0
## CD28 0 2.038401 0.553 0.012 0
## CTLA4 0 1.319968 0.289 0.001 0
## ICOS 0 2.520687 0.539 0.002 0
## IKZF2 0 2.335888 0.382 0.001 0
## TNIK 0 3.303764 0.632 0.002 0
## GPRIN3 0 2.874834 0.632 0.014 0
## LEF1 0 1.911537 0.263 0.001 0
## INPP4B 0 3.156366 0.671 0.002 0
## IL7R 0 1.925463 0.395 0.001 0
## GZMK 0 1.400449 0.197 0.000 0
## PAM 0 2.577763 0.645 0.002 0
## CAMK4 0 2.665704 0.487 0.002 0
## ITK 0 2.808106 0.645 0.008 0
## AC010609.1 0 1.117180 0.250 0.001 0
## THEMIS 0 2.183010 0.289 0.001 0
## TRBC1 0 1.460282 0.224 0.001 0
With all the evidence gathered, we proceed to exclude cluster 1 and 8:
se_l$M102 <- se_l$M102[, !(se_l$M102$seurat_clusters %in% c("1", "8"))]
DimPlot(se_l$M102, cols = cols)
We now reprocess the dataset and cluster cells into meaningful cell states:
se_l$M102 <- process_seurat(se_l$M102)
DimPlot(se_l$M102, cols = cols)
We know from exploratory analysis that the main driver of intratumoral transcriptional heterogeneity in MCL are copy number alterations. In particular, we detected a loss of chromosome Y (among others) using infercnv. We visualize the expression of genes encoded in the Y chromosome, collapse their expression into a per cell score (“chr Y score”) and classify them in chr Y+ or chr Y- using this score:
# Visualize expression of genes encoded in chromosome Y
goi_chrY <- c("UTY", "KDM5D", "DDX3Y", "USP9Y", "ZFY", "EIF1AY")
FeaturePlot(se_l$M102, features = goi_chrY) &
scale_color_viridis_c(option = "magma") &
NoLegend()
# Score cells based on chromosome Y expression
se_l$M102 <- AddModuleScore(
se_l$M102,
features = list(goi_chrY),
name = "chrY_score"
)
FeaturePlot(se_l$M102, "chrY_score1") &
scale_color_viridis_c(option = "magma")
# Classify cells in chr Y + and chr Y -
(density_gg <- se_l$M102@meta.data %>%
ggplot(aes(chrY_score1)) +
geom_density() +
geom_vline(xintercept = chrY_cutoff_102, linetype = "dashed", color = "darkblue") +
theme_classic())
se_l$M102$has_loss_chrY <- ifelse(
se_l$M102$chrY_score1 > chrY_cutoff_102,
"chrY+",
"chrY-"
)
DimPlot(se_l$M102, group.by = "has_loss_chrY")
Strategy: identify broad clusters (chrY+ and -), subcluster inside each of them with FindSubCluster:
se_l$M102 <- FindClusters(se_l$M102, resolution = 0.05)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8499
## Number of edges: 308665
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9663
## Number of communities: 3
## Elapsed time: 0 seconds
DimPlot(se_l$M102)
# Subcluster cluster 0 (chr Y+)
se_l$M102 <- FindSubCluster(
se_l$M102,
graph.name = "RNA_snn",
resolution = 0.6,
cluster = "0",
subcluster.name = "subclusters_0"
)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 5901
## Number of edges: 209385
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7618
## Number of communities: 7
## Elapsed time: 0 seconds
DimPlot(se_l$M102, group.by = "subclusters_0")
# Subcluster cluster 1 (chr Y-)
se_l$M102$clusters_20230625 <- se_l$M102$subclusters_0
Idents(se_l$M102) <- "clusters_20230625"
se_l$M102 <- FindSubCluster(
se_l$M102,
graph.name = "RNA_snn",
resolution = 0.6,
cluster = "1",
subcluster.name = "clusters_20230625"
)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2292
## Number of edges: 85715
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7027
## Number of communities: 5
## Elapsed time: 0 seconds
DimPlot(se_l$M102, group.by = "clusters_20230625")
# Further subcluster to get non-tumoral cells
Idents(se_l$M102) <- "clusters_20230625"
se_l$M102 <- FindSubCluster(
se_l$M102,
graph.name = "RNA_snn",
resolution = 0.1,
cluster = "1_4",
subcluster.name = "clusters_20230625"
)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 80
## Number of edges: 1431
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9126
## Number of communities: 2
## Elapsed time: 0 seconds
DimPlot(se_l$M102, group.by = "clusters_20230625")
Idents(se_l$M102) <- "clusters_20230625"
# Markers subclusters cluster 0
se_clus0 <- se_l$M102[, str_detect(se_l$M102$clusters_20230625, "^0")]
DimPlot(se_clus0)
markers_0 <- FindAllMarkers(se_clus0, only.pos = TRUE, logfc.threshold = 0.75)
markers_0 <- markers_0[markers_0$p_val_adj < 0.001, ]
DT::datatable(markers_0, options = list(scrollX = TRUE))
# Markers subclusters cluster 1
se_clus1 <- se_l$M102[, str_detect(se_l$M102$clusters_20230625, "^1")]
Idents(se_clus1) <- "clusters_20230625"
DimPlot(se_clus1)
markers_1 <- FindAllMarkers(se_clus1, only.pos = TRUE, logfc.threshold = 0.75)
markers_1 <- markers_1[markers_1$p_val_adj < 0.001, ]
DT::datatable(markers_1, options = list(scrollX = TRUE))
Idents(se_l$M102) <- "clusters_20230625"
Annotation:
| Subcluster | Markers | Annotation |
|---|---|---|
| 0_0 | TCF4, RGS12 | c1 TCF4Hi RGS12Hi |
| 0_1 | MT1G, MT1E, MT1X, MT2A | c1 metallothionein |
| 0_2 | MT1G, MT1E, MT1X, MT2A | c1 metallothionein |
| 0_3 | MYC, NFKB1, MIR155HG | c1 MYC+NFKB1+MIR155HG+ |
| 0_4 | CD69, JUN, FOS | c1 CD69+AP1+ |
| 0_5 | TUBB, CENPP | c1 S phase |
| 0_6 | CCL22, NME1 | c1 CCL22+NME1+ |
| 1_0 | MT1G, MT1E, MT1X, MT2A | c2 metallothionein |
| 1_1 | CD69, JUN, FOS | c2 CD69+AP1+ |
| 1_2 | TCF4, RGS12 | c2 TCF4Hi RGS12Hi |
| 1_3 | MYC, NFKB1, MIR155HG | c2 MYC+NFKB1+MIR155HG+ |
| 1_4_0 | NR3C2 | c2 NR3C2+ |
| 1_4_1 | SOX11-, CCND1- | c3 non-tumoral B-cells |
| 2 | MKI67, CENPF | c1 G2M Phase |
current_levels <- c("0_0", "0_1", "0_2", "0_3", "0_4", "0_5", "0_6", "1_0", "1_1",
"1_2", "1_3", "1_4_0", "1_4_1", "2")
new_levels <- c("c1 TCF4Hi RGS12Hi", "c1 metallothionein", "c1 metallothionein",
"c1 MYC+NFKB1+MIR155HG+", "c1 CD69+AP1+", "c1 S phase",
"c1 CCL22+NME1+", "c2 metallothionein", "c2 CD69+AP1+",
"c2 TCF4Hi RGS12Hi", "c2 MYC+NFKB1+MIR155HG+", "c2 NR3C2+",
"c3 non-tumoral B-cells", "c1 G2M Phase")
names(new_levels) <- current_levels
se_l$M102$annotation_20230625 <- new_levels[se_l$M102$clusters_20230625]
Idents(se_l$M102) <- "annotation_20230625"
cols <- kelly.colors(length(new_levels) + 1)
names(cols) <- NULL
cols <- cols[2:length(cols)]
DimPlot(se_l$M102, pt.size = 0.85, cols = cols)
Plot markers
goi <- c(goi_chrY, "TCF4", "RGS12", "MT1G", "MT1E", "MT1X", "MT2A", "MYC",
"NFKB1", "MIR155HG", "CD69", "JUN", "FOS", "TUBB", "CENPP", "CCL22",
"NME1", "NR3C2", "SOX11", "CCND1", "MKI67", "CENPF")
dot_plot <- DotPlot(se_l$M102, features = rev(goi)) +
scale_color_distiller(palette = "Blues", direction = 1) +
theme(
axis.title = element_blank(),
axis.text.y = element_text(size = 8),
legend.text = element_text(size = 7),
legend.title = element_text(size = 7),
axis.text.x = element_text(size = 8, angle = 90, hjust = 1, vjust = 1)
)
dot_plot
We follow the same process for the second donor (M413):
se_l$M413 <- FindClusters(se_l$M413, resolution = 0.4)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11074
## Number of edges: 399988
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8898
## Number of communities: 13
## Elapsed time: 1 seconds
cols <- kelly.colors(length(levels(se_l$M413$seurat_clusters)) + 1)
names(cols) <- NULL
cols <- cols[2:length(cols)]
vln_qc_clusters <- VlnPlot(
se_l$M413,
c("nCount_RNA", "nFeature_RNA", "pct_mt", "doublet_score"),
pt.size = 0,
ncol = 2,
cols = cols
)
vln_qc_clusters[[1]] <- vln_qc_clusters[[1]] +
scale_y_log10()
vln_qc_clusters
umap_dbl_score <- FeaturePlot(se_l$M413, "doublet_score") +
scale_color_viridis_c(option = "magma")
umap_clusters <- DimPlot(se_l$M413, group.by = "seurat_clusters", cols = cols)
umap_clusters | umap_dbl_score # cluster 4 seems to be composed of doublets
FeaturePlot(se_l$M413, c("CD3E", "CD3D")) # cluster 9 might be T cells
markers_4_9 <- map(c("4", "9"), \(x) {
df <- FindMarkers(
se_l$M413,
ident.1 = x,
only.pos = TRUE
)
head(df, 40)
})
map(markers_4_9, head, 30)
## [[1]]
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## NEB 1.489423e-193 0.8185946 0.885 0.335 5.451437e-189
## RASSF6 6.840967e-192 0.7913718 0.662 0.190 2.503862e-187
## LINC01170 3.009270e-187 0.7972957 0.446 0.097 1.101423e-182
## PCDH9 4.679998e-184 0.5076865 0.594 0.156 1.712926e-179
## CLNK 2.964332e-169 0.5101595 0.615 0.173 1.084975e-164
## LINC01876 2.604968e-168 0.4905777 0.426 0.094 9.534442e-164
## ZNF331 9.665181e-154 0.9075758 0.669 0.242 3.537553e-149
## SAMD12 7.070175e-150 0.4835986 0.706 0.234 2.587755e-145
## LHFPL2 9.661333e-146 0.6648814 0.952 0.460 3.536145e-141
## AC011029.1 2.085981e-138 0.4222857 0.275 0.050 7.634900e-134
## CHPT1 3.015170e-129 0.7396944 0.740 0.288 1.103582e-124
## LINC01811 6.944307e-129 0.4236081 0.425 0.114 2.541686e-124
## DIP2C 4.072152e-124 0.4898166 0.622 0.214 1.490448e-119
## IFNG-AS1 1.991248e-123 0.5814476 0.925 0.454 7.288167e-119
## AL589693.1 1.471374e-120 0.4877610 0.702 0.266 5.385377e-116
## POF1B 7.695993e-119 0.2904005 0.479 0.141 2.816810e-114
## KHDRBS3 6.644928e-117 0.4994470 0.250 0.049 2.432110e-112
## CRY1 8.937637e-114 0.3976285 0.455 0.137 3.271265e-109
## ZNF10 2.886105e-113 0.5467851 0.613 0.224 1.056343e-108
## LINC01949 9.888942e-109 0.2903443 0.291 0.067 3.619452e-104
## NR4A2 7.925605e-108 0.4388300 0.410 0.122 2.900851e-103
## CHL1 2.644286e-106 0.2986969 0.307 0.075 9.678352e-102
## BCL2L11 4.110751e-104 0.4535714 0.474 0.156 1.504576e-99
## PPP1R9A 4.373997e-104 0.3943337 0.288 0.068 1.600927e-99
## IRS2 7.435251e-102 0.2777543 0.208 0.039 2.721376e-97
## RGS2 1.814152e-101 0.4177176 0.446 0.145 6.639977e-97
## TNFRSF21 6.634843e-101 0.2778408 0.232 0.048 2.428419e-96
## C12orf74 1.370321e-99 0.3635609 0.549 0.193 5.015511e-95
## NR4A3 2.671434e-99 0.2874947 0.309 0.080 9.777714e-95
## KCNQ1 1.571569e-97 0.2772179 0.468 0.152 5.752099e-93
##
## [[2]]
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## CD2 0 1.6735759 0.419 0.004 0
## CD247 0 2.5755216 0.648 0.011 0
## HNRNPLL 0 1.9191762 0.512 0.011 0
## AAK1 0 2.2796268 0.615 0.038 0
## LINC01934 0 3.2191254 0.548 0.041 0
## PLCL1 0 2.9147021 0.505 0.017 0
## CD28 0 1.1446739 0.289 0.002 0
## ICOS 0 1.4612454 0.316 0.003 0
## TRAT1 0 0.8116984 0.186 0.002 0
## CD96 0 2.5146420 0.714 0.029 0
## TNIK 0 2.3195393 0.538 0.020 0
## DTHD1 0 1.7384070 0.329 0.004 0
## PTPN13 0 1.3406290 0.226 0.004 0
## GPRIN3 0 2.3136307 0.575 0.017 0
## LEF1 0 1.6233189 0.302 0.004 0
## INPP4B 0 3.2944994 0.718 0.013 0
## AC139720.1 0 1.6323853 0.229 0.004 0
## IL7R 0 1.1029533 0.226 0.002 0
## FYB1 0 3.1914259 0.827 0.039 0
## PAM 0 2.4210930 0.585 0.027 0
## CAMK4 0 1.8793804 0.405 0.015 0
## PPP2R2B 0 1.8381960 0.299 0.006 0
## ITK 0 1.8560248 0.482 0.010 0
## LCP2 0 1.8670422 0.548 0.025 0
## THEMIS 0 3.4031080 0.721 0.011 0
## TOX 0 3.1708965 0.801 0.083 0
## TTC39B 0 1.6936300 0.419 0.017 0
## MLLT3 0 2.1506749 0.472 0.031 0
## PRKCQ 0 1.5945736 0.415 0.006 0
## ST8SIA1 0 3.3122906 0.751 0.011 0
# scDblFinder
sce_M413 <- as.SingleCellExperiment(se_l$M413)
dbl_out <- findDoubletClusters(
sce_M413,
clusters = sce_M413$seurat_clusters
)
dbl_out
## DataFrame with 13 rows and 9 columns
## source1 source2 num.de median.de best p.value lib.size1 lib.size2 prop
## <character> <character> <integer> <numeric> <character> <numeric> <numeric> <numeric> <numeric>
## 0 12 5 13 463.0 HSP90AA1 1.64605e-09 1.009592 1.785401 0.31497201
## 12 11 3 40 303.5 HSPA1B 1.26195e-33 1.820129 1.063446 0.00758534
## 2 12 11 46 603.5 TSHZ2 1.06602e-41 1.002583 1.824830 0.10032509
## 3 12 4 67 974.5 HSPA1B 5.28304e-16 0.940339 2.564506 0.08994040
## 11 12 2 100 381.5 SKAP1 1.28777e-21 0.549412 0.547996 0.01354524
## ... ... ... ... ... ... ... ... ... ...
## 4 12 1 409 1436.0 CD58 1.60262e-09 0.366675 0.418753 0.0831678
## 10 12 9 413 1184.5 FNDC3B 6.05853e-53 0.272392 0.321256 0.0135452
## 6 12 8 440 1359.5 MIR155HG 1.27779e-35 0.773832 1.501558 0.0469568
## 9 11 10 540 1126.5 PRKCH 6.34361e-90 1.543282 3.112780 0.0271808
## 8 12 7 1357 2414.5 AC023590.1 9.99730e-146 0.515353 0.545560 0.0273614
# Exclude cluster 4 (doublets) and 9 (T cells) and reprocess
se_l$M413 <- se_l$M413[, !(se_l$M413$seurat_clusters %in% c("4", "9"))]
DimPlot(se_l$M413, cols = cols)
se_l$M413 <- process_seurat(se_l$M413)
DimPlot(se_l$M413, cols = cols)
Similar as before, we classify cells into chrY+ and chrY-:
# Visualize expression of genes encoded in chromosome Y
FeaturePlot(se_l$M413, features = goi_chrY) &
scale_color_viridis_c(option = "magma") &
NoLegend()
# Score cells based on chromosome Y expression
se_l$M413 <- AddModuleScore(
se_l$M413,
features = list(goi_chrY),
name = "chrY_score"
)
FeaturePlot(se_l$M413, "chrY_score1") &
scale_color_viridis_c(option = "magma")
# Classify cells in chr Y + and chr Y -
(density_gg <- se_l$M413@meta.data %>%
ggplot(aes(chrY_score1)) +
geom_density() +
geom_vline(xintercept = chrY_cutoff_413, linetype = "dashed", color = "darkblue") +
theme_classic())
se_l$M413$has_loss_chrY <- ifelse(
se_l$M413$chrY_score1 > chrY_cutoff_413,
"chrY+",
"chrY-"
)
DimPlot(se_l$M413, group.by = "has_loss_chrY")
Cluster:
se_l$M413 <- FindClusters(se_l$M413, resolution = 0.05)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 9852
## Number of edges: 366652
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9702
## Number of communities: 5
## Elapsed time: 1 seconds
DimPlot(se_l$M413)
# Subcluster cluster 0 (chr Y+)
se_l$M413 <- FindSubCluster(
se_l$M413,
graph.name = "RNA_snn",
resolution = 0.5,
cluster = "0",
subcluster.name = "subclusters_0"
)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 6709
## Number of edges: 241651
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8174
## Number of communities: 8
## Elapsed time: 0 seconds
DimPlot(se_l$M413, group.by = "subclusters_0")
# Further subcluster cluster 0_2
Idents(se_l$M413) <- "subclusters_0"
se_l$M413 <- FindSubCluster(
se_l$M413,
graph.name = "RNA_snn",
resolution = 0.2,
cluster = "0_2",
subcluster.name = "clusters_20230625"
)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1103
## Number of edges: 35158
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8297
## Number of communities: 2
## Elapsed time: 0 seconds
DimPlot(se_l$M413, group.by = "clusters_20230625")
Idents(se_l$M413) <- "clusters_20230625"
# Finally, let us further subcluster cluster 3, which contains normal GCBC
se_l$M413 <- FindSubCluster(
se_l$M413,
graph.name = "RNA_snn",
resolution = 0.1,
cluster = "2",
subcluster.name = "clusters_20230625"
)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 556
## Number of edges: 17577
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9406
## Number of communities: 2
## Elapsed time: 0 seconds
cols <- kelly.colors(length(unique(se_l$M413$clusters_20230625)) + 1)
names(cols) <- NULL
cols <- cols[2:length(cols)]
DimPlot(se_l$M413, group.by = "clusters_20230625", cols = cols)
Idents(se_l$M413) <- "clusters_20230625"
markers_413 <- FindAllMarkers(se_l$M413, only.pos = TRUE, logfc.threshold = 0.75)
markers_413 <- markers_413[markers_413$p_val_adj < 0.001, ]
DT::datatable(markers_413, options = list(scrollX = TRUE))
Annotation:
| Subcluster | Markers | Annotation |
|---|---|---|
| 0_0 | chrY-, CD69lo, CXCR4lo, JUNB- | c1 chrY-CD69LoCXCR4LoJUNB- |
| 0_1 | chrY+, CD69lo, CXCR4lo, JUNB- | c1 chrY+CD69LoCXCR4LoJUNB- |
| 0_2_0 | CD69, CXCR4, JUNB | c1 CD69Hi CXCR4HiJUNB+ |
| 0_2_1 | CD69, CXCR4, JUNB, TSHZ2, MECOM | c1 CD69Hi CXCR4HiJUNB+TSHZ2+MECOM+ |
| 0_3 | TSHZ2, MECOM | c1 TSHZ2+MECOM+ |
| 0_4 | CENPK, POLQ | c1 S phase |
| 0_5 | MYC, MIR155HG, NFKB1 | c1 MYC+MIR155HG+NFKB1Hi |
| 0_6 | MT-ND2, MT-ND3, MT-ATP8 | poor-quality |
| 0_7 | HSPD1, HSP90AB1, HSPA1A | poor-quality |
| 1 | PCDH9 | c2 PCDH9 Hi |
| 2_0 | MEF2B, RGS13 | GCBC |
| 2_1 | CENPF, CENPE | c1 G2M phase |
| 3 | CD96, TOX, SOX5 | NBC/MBC |
| 4 | XBP1, PRDM1 | PC |
current_levels <- c("0_0", "0_1", "0_2_0", "0_2_1", "0_3", "0_4", "0_5", "0_6",
"0_7", "1", "2_0", "2_1", "3", "4")
new_levels <- c("c1 chrY-CD69LoCXCR4LoJUNB-", "c1 chrY+CD69LoCXCR4LoJUNB-",
"c1 CD69Hi CXCR4HiJUNB+", "c1 CD69Hi CXCR4HiJUNB+TSHZ2+MECOM+",
"c1 TSHZ2+MECOM+", "c1 S phase", "c1 MYC+MIR155HG+NFKB1Hi",
"poor-quality", "poor-quality", "c2 PCDH9 Hi", "GCBC", "c1 G2M phase",
"NBC/MBC", "PC")
names(new_levels) <- current_levels
se_l$M413$annotation_20230625 <- new_levels[se_l$M413$clusters_20230625]
Idents(se_l$M413) <- "annotation_20230625"
cols <- kelly.colors(length(new_levels) + 1)
names(cols) <- NULL
cols <- cols[2:length(cols)]
DimPlot(se_l$M413, pt.size = 0.85, cols = cols)
# Remove cluster 6 and reprocess
se_l$M413 <- se_l$M413[, se_l$M413$annotation_20230625 != "poor-quality"]
se_l$M413 <- process_seurat(se_l$M413)
DimPlot(se_l$M413, pt.size = 0.85, cols = cols)
Plot markers
goi_413 <- c("CCND1", "SOX11", "CD69", "FOS", "ZNF331", "RGS2" , "MECOM",
"TSHZ2", "JUNB", "CXCR4", "CCL3", "PCDH9", "MIR155HG", "NFKB1",
"PCNA", "TOP2A", "MCM6", "MKI67", "MYC", "IRF4", "PRDM1", "BCL6",
"MEF2B", "PLAC8", "SOX5")
dot_plot <- DotPlot(se_l$M413, features = rev(goi_413)) +
scale_color_distiller(palette = "Blues", direction = 1) +
theme(
axis.title = element_blank(),
axis.text.y = element_text(size = 8),
legend.text = element_text(size = 7),
legend.title = element_text(size = 7),
axis.text.x = element_text(size = 8, angle = 90, hjust = 1, vjust = 1)
)
dot_plot
dir.create(here("MCL/results/R_objects"), recursive = TRUE)
saveRDS(se_l, here("MCL/results/R_objects/20230625_seurat_list_MCL_M102_M413_tumoral_cells_annotated.rds"))
sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /home/groups/singlecell/rmassoni/anaconda3/envs/richter/lib/libopenblasp-r0.3.21.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Polychrome_1.5.1 glue_1.6.2 here_1.0.1 ggpubr_0.6.0 DT_0.27 scDblFinder_1.12.0 SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0 Biobase_2.58.0 GenomicRanges_1.50.0 GenomeInfoDb_1.34.9 IRanges_2.32.0 S4Vectors_0.36.0 BiocGenerics_0.44.0 MatrixGenerics_1.10.0 matrixStats_0.63.0 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.0 purrr_1.0.1 readr_2.1.4 tidyr_1.3.0 tibble_3.1.8 ggplot2_3.4.1 tidyverse_1.3.2 SeuratObject_4.1.3 Seurat_4.3.0 BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] rtracklayer_1.58.0 scattermore_0.8 R.methodsS3_1.8.2 bit64_4.0.5 knitr_1.42 R.utils_2.12.2 irlba_2.3.5.1 DelayedArray_0.24.0 data.table_1.14.6 RCurl_1.98-1.10 generics_0.1.3 ScaledMatrix_1.6.0 cowplot_1.1.1 RANN_2.6.1 future_1.31.0 bit_4.0.5 tzdb_0.3.0 spatstat.data_3.0-0 xml2_1.3.3 lubridate_1.9.1 httpuv_1.6.9 assertthat_0.2.1 viridis_0.6.2 gargle_1.3.0 xfun_0.37 hms_1.1.2 jquerylib_0.1.4 evaluate_0.20 promises_1.2.0.1 fansi_1.0.4 restfulr_0.0.15 dbplyr_2.3.0 readxl_1.4.2 igraph_1.3.5 DBI_1.1.3 htmlwidgets_1.6.1 spatstat.geom_3.0-6 googledrive_2.0.0 ellipsis_0.3.2 crosstalk_1.2.0 backports_1.4.1 bookdown_0.33 deldir_1.0-6 sparseMatrixStats_1.10.0 vctrs_0.5.2
## [46] ROCR_1.0-11 abind_1.4-5 cachem_1.0.6 withr_2.5.0 progressr_0.13.0 vroom_1.6.1 sctransform_0.3.5 GenomicAlignments_1.34.0 scran_1.26.2 goftest_1.2-3 cluster_2.1.4 lazyeval_0.2.2 crayon_1.5.2 spatstat.explore_3.0-6 labeling_0.4.2 edgeR_3.40.2 pkgconfig_2.0.3 nlme_3.1-162 vipor_0.4.5 rlang_1.0.6 globals_0.16.2 lifecycle_1.0.3 miniUI_0.1.1.1 modelr_0.1.10 rsvd_1.0.5 ggrastr_1.0.1 cellranger_1.1.0 rprojroot_2.0.3 polyclip_1.10-4 lmtest_0.9-40 Matrix_1.5-3 carData_3.0-5 zoo_1.8-11 reprex_2.0.2 beeswarm_0.4.0 ggridges_0.5.4 googlesheets4_1.0.1 png_0.1-8 viridisLite_0.4.1 rjson_0.2.21 bitops_1.0-7 R.oo_1.25.0 KernSmooth_2.23-20 Biostrings_2.66.0 DelayedMatrixStats_1.20.0
## [91] parallelly_1.34.0 spatstat.random_3.1-3 rstatix_0.7.2 ggsignif_0.6.4 beachmat_2.14.2 scales_1.2.1 magrittr_2.0.3 plyr_1.8.8 ica_1.0-3 zlibbioc_1.44.0 compiler_4.2.2 dqrng_0.3.0 BiocIO_1.8.0 RColorBrewer_1.1-3 fitdistrplus_1.1-8 Rsamtools_2.14.0 cli_3.6.0 XVector_0.38.0 listenv_0.9.0 patchwork_1.1.2 pbapply_1.7-0 MASS_7.3-58.2 tidyselect_1.2.0 stringi_1.7.12 highr_0.10 yaml_2.3.7 BiocSingular_1.14.0 locfit_1.5-9.7 ggrepel_0.9.3 grid_4.2.2 sass_0.4.5 tools_4.2.2 timechange_0.2.0 future.apply_1.10.0 parallel_4.2.2 bluster_1.8.0 metapod_1.6.0 gridExtra_2.3 farver_2.1.1 scatterplot3d_0.3-44 Rtsne_0.16 digest_0.6.31 BiocManager_1.30.19 shiny_1.7.4 Rcpp_1.0.10
## [136] car_3.1-2 broom_1.0.3 scuttle_1.8.4 later_1.3.0 RcppAnnoy_0.0.20 httr_1.4.4 colorspace_2.1-0 rvest_1.0.3 XML_3.99-0.13 fs_1.6.1 tensor_1.5 reticulate_1.26 splines_4.2.2 uwot_0.1.14 statmod_1.5.0 spatstat.utils_3.0-1 scater_1.26.1 sp_1.6-0 xgboost_1.7.5.1 plotly_4.10.1 xtable_1.8-4 jsonlite_1.8.4 R6_2.5.1 pillar_1.8.1 htmltools_0.5.4 mime_0.12 fastmap_1.1.0 BiocParallel_1.32.5 BiocNeighbors_1.16.0 codetools_0.2-19 utf8_1.2.3 lattice_0.20-45 bslib_0.4.2 spatstat.sparse_3.0-0 ggbeeswarm_0.7.1 leiden_0.4.3 survival_3.5-3 limma_3.54.0 rmarkdown_2.20 munsell_0.5.0 GenomeInfoDbData_1.2.9 haven_2.5.1 reshape2_1.4.4 gtable_0.3.1